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Summary 

A numerical method for the convective heat transfer 
problem is developed for low speed flow at mild 
temperatures. A simplified energy equation is added to 
the incompressible Navier-Stokes formulation by using 
Boussinesq approximation to account for the buoyancy 
force. A pseudocompressibility method is used to solve 
the resulting set of equations for steady-state solutions in 
conjunction with an approximate factorization scheme. 

A Neumann-type pressure boundary condition is devised 
to account for the interaction between pressure and 
temperature terms, especially near a heated or cooled 
solid boundary. It is shown that the present method is 
capable of predicting the temperature field in an 
incompressible flow. 

Introduction 

Heat transfer in viscous incompressible flow is of interest 
in many industrial applications. For example, in a liquid 
rocket engine, the liquid fuel and oxidizer are used as 
coolant in various components such as the bearing in 
the oxidizer and fuel turbopump. Also, the flow in an 
autoclave for curing aerospace parts can be analyzed 
using an incompressible flow assumption. For a complete 
analysis of heat transfer in a wide range of temperature, 
one must include radiation effects as well as boiling heat 
transfer. However, of current interest are the problems 
dominated by convective heat transfer. Therefore, in the 
present study, the internal energy generated by viscous 
dissipation and the thermal radiation effects are neglected. 
The fluid is assumed to be incompressible with constant 
physical properties except for the buoyancy effect due to 
density variations. When the temperature of the flow field 
is not high, the thermally driven velocity is small relative 
to sonic speed. Thus a Boussinesq approximation can be 
applied to the incompressible Navier-Stokes equations 
to represent the temperature field. The purpose of the 
present study is to develop a computational capability for 
simulating viscous incompressible flows with temperature 
variations. Since the method is intended for application 
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to three-dimensional problems, a primitive variable 
formulation is chosen based on a structured-grid 
approach. To use one of the primitive variable solvers 
(refs. 1-5), a simplified energy equation is added to the 
incompressible Navier-Stokes equations. In the present 
study, the first incompressible Navier-Stokes solver 
developed at Ames Research Center, the INS3D code 
(ref. 1), is selected to test the feasibility of using the 
present formulation in predicting the temperature field 
in an incompressible medium. 

To validate the flow solver, a simple channel flow is 
computed where an analytical solution exists. Then, two- 
dimensional flow problems are computed and compared 
with other numerical and experimental results. In all 
these problems, natural, mixed, and forced convection 
problems are examined. Finally, computed results in 
three dimensions are compared with experiments. The 
simulation capability related to thermal effects has been 
demonstrated. 

Solution Methods 

Boussinesq Approximation 

Neglecting the adiabatic temperature increase due 
to friction, the equations governing the flow of an 
incompressible fluid with constant properties can be 
written as 
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where xj is the Cartesian coordinates, Uj the correspond- 
ing velocity components, p the pressure, t the time, Tjj 
the viscous-stress tensor, g the vector of gravitational 
acceleration, p the density, T the temperature, and a the 
thermal diffusivity. The viscous stress tensor can be 
written as 
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where v is the kinematic viscosity, Sy the strain-rate 
tensor, and Ry the Reynolds stresses. Various levels of 
closure models for Ry are possible. In the present study, 
turbulence is simulated by an eddy viscosity model using 
a constitutive equation of the following form: 

Rij=jRkk8ij-2v t Sij (6) 

where v t is the turbulent eddy viscosity. By including the 
normal stress, R]±, in the pressure, v in equation (4) can 
be replaced by (v + v t ) as follows: 

^ij = 2(v + v t )Sy = 2v t Sy (7) 


In the remainder of this report, the total viscosity, vj, will 
be represented simply by v. The present formulations 
allow for a spatially varying viscosity. 

The buoyancy force term is simplified through the 
Boussinesq approximation where the density in the 
buoyancy term is represented by a linear variation of the 
temperature, 

p = Po{l-p(T-T 0 )} (8) 


where P is the coefficient of thermal expansion. The 
buoyancy term based on this approximation is included in 
the momentum equation (2) resulting in 
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The above governing equations can be nondimensional- 
ized by introducing the following dimensionless 
quantities: 
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Here, uq is the reference velocity, L the reference length, 
and Ti - To the reference temperature difference. By 
omitting the prime in the nondimensional variables, the 
governing equations with Boussinesq approximation can 
be written as the following dimensionless form: 
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Where, e„ is the unit vector for gravitational accelera- 
tion. Depending on the flow regime, the reference 
quantities vary and the coefficients are defined 
accordingly: 

for natural convection. 


Cm = — » C B = — , Cp = — — ■ 
M Re Re 2 RePr 

for forced/mixed convection, 

Cm = Pr , C B = Ra Pr , Cg = 1 
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Here, the nondimensional numbers are defined as 
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Richardson number: 



for natural convection, 
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for forced/mixed convection, 
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Pseudocompressibility Formulation in Generalized 
Coordinates 

The pseudocompressibility is introduced after the 
governing equations (9)-(l 1) are transformed into general 
curvilinear coordinates, (^,r|,Q, which results in 
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J = Jacobian of the transformation 

(3 = pseudocompressibility parameter 

Here, X\ , A, 2 , and X 3 are components of the unit tensor e g 
in the x, y, and z directions, respectively. 
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H = Jacobian matrix of the source term S 
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(r 2 and T 3 are similarly defined) 


Sp = finite difference form of — 
a = 1/2 for trapezoidal, or 1 for Euler 


l m =diag[0, 1,1,1] 


This equation is iterated in pseudotime until the solution 
converges to steady state, at which time the original 
incompressible Navier-Stokes equations are satisfied. 

A direct inversion of equation (15) would become a 
Newton iteration for a steady-state solution. In three 
dimensions, however, direct inversion of a large block 
banded matrix of the unfactored scheme would be 
impractical. Numerous iterative schemes can be imple- 
mented to solve these equations (see ref. 6 for a review). 
In the present study, an approximate factorization scheme 
by Beam and Warming (ref. 7) is used. 


Numerical Method 


Buoyancy Effect on Pressure 


An unfactored implicit scheme can be obtained by 
linearizing the flux vectors with respect to the previous 
time step and dropping terms of the second and higher 
order, which results in the following equations in delta 
form: 
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The buoyancy effect in a thermally convective flow needs 
to be assessed relative to the pressure wave propagation 
and the boundary layer development. For the pseudo- 
compressibility formulation, the pressure wave propa- 
gates at a finite speed, the magnitude of which depends 
on the pseudocompressibility parameter. When the 
thermal effect is the dominant driving force such as in 
natural convection, a pressure gradient is created by the 
temperature variations. Thus the pressure boundary 
condition should include temperature effect. A full 
account of buoyancy effect on pseudocompressibility 
will be given in a later report. 


3 



Computed Results 

A Vertical Channel with Temperature Gradient 


A two-dimensional vertical channel flow is considered as 
a first test case of the present formulation. As shown in 
figure 1(a), one wall of the channel is heated to T w and 
the other is cooled down to -T w . For an infinitely long 
channel, an analytic solution exists in the following form: 


v = 6V 0 ifl 





T = T »( 1 - 2 f) 


where v is the vertical velocity, h the channel width, and 
Vo an average velocity. 

The computation was performed using a long channel 
with L = lOOh. At the inlet, a uniform temperature T = 0 
and a constant velocity v =Vo are specified, while at the 
outlet a Neumann condition is imposed. In a forced 
convection mode, the flow becomes a two-dimensional 
Poiseuille flow. In the case of natural convection, which 
has the zero streamwise pressure gradient, the flow is 
generated by the buoyancy force, whereas in mixed 
convection the flow develops not only by the buoyancy 
force but also by the streamwise pressure gradient. All 
three modes of heat transfer problems are computed by 
selecting nondimensional parameters to represent 
respective flow fields. In figures l(b)-l(d), computed 
temperature, vorticity, and velocity profiles at y/h = 50 
are compared with the analytic solutions. Computations 
essentially reproduced the analytic solutions. 


Flow around a Heated Circular Cylinder 

The external flow test case consisted of a circular cylinder 
under the influence of a uniform upwardly moving fluid. 
In two dimensions, the stream function-vorticity formu- 
lation has been used in numerous numerical studies (for 
example, see ref. 8). Since detailed measurements of the 
flow field involving heat transfer are rare, the results of 
the present incompressible Navier-Stokes computation 
are compared to those of stream function-vorticity 
approach by Sa (ref. 8). In the figure, INS and SV 
represent the results obtained using incompressible 
Navier-Stokes and stream function-vorticity formula- 
tions, respectively. Since the stream function-vorticity 
formulation satisfies the divergence free velocity 
condition, this comparison can also be used to show that 
the present pseudocompressible formulation results in 
incompressible solutions. 


In figures 2 and 3, computed results for forced, mixed, 
and natural convection cases are presented. The incom- 
pressible Navier-Stokes results are shown on the left half 
of each figure and results of the stream function-vorticity 
are shown on the right half (see fig. 2). The forced con- 
vection case was computed at Reynolds numbers ranging 
from 5 to 40. However, in the present paper, only the 
results for Re = 20 are reported. The mixed convection 
for the heated or cooled cylinder was computed at the 
Richardson number, which is defined in equation (13), 
from -1 to 4 with Re = 20. When the Richardson number 
is negative, the cylinder is cooled, and the fluid in the 
boundary layer and in the wake region is decelerated by 
the cooling. As the Richardson number increases, the 
flow is accelerated and the separation of the boundary 
layer is suppressed. The streamlines in figure 2(b) 
indicate that there is no separation at Ri = 4. The natural 
convection case was computed at the Rayleigh number up 
to 105. The dimensionless heat transfer coefficient, the 
Nusselt number, is compared with experiments and other 
computations as shown in figure 3 for the forced, mixed, 
and natural convection cases. As the Reynolds number 
increases or the cylinder is heated, the stronger velocity 
near a surface makes the Nusselt number increase. 
Overall, the present results agree well with numerical 
and experimental data reported in references 9-14. 


Thermally Driven Bifurcation in a Rectangular Cavity 

Thermal instability is investigated next for a rectangular 
cavity with an aspect ratio of 1 or 2, which has a hot 
bottom wall and a cold top wall. For an aspect ratio of 1 .0 
at Ra = 105 , there exists a unique solution with a single 
vortex as shown in figure 4(a). The INS (left) and the SV 
(right) show excellent agreement. On the other hand, the 
cavity with an aspect ratio of 2.0 at Ra = 10^ may have 
two types of solutions: a double vortex as shown in 
figure 4(b) or a single vortex in the middle of the cavity 
as shown in figure 4(c). The bifurcation depends on the 
external disturbances and the initial and boundary 
conditions. In figure 4(b), the result of the INS (left) is a 
little different from that of the SV (right), since the grid is 
too coarse near the center line (the same grid number of 
21x21 was used for both aspect ratios). However, in 
general it is shown that the present method is adequate to 
simulate the thermal instability. 


Three-Dimensional Cavity 

Most three-dimensional experiments are focused on the 
investigation of heat transfer coefficients at surfaces for 
practical applications. However, Morrison and Tran 
(ref. 14) experimentally investigated the flow structure 
in a natural convection mode generated by heated walls 
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in a vertical rectangular cavity shown in figure 5. The 
temperature difference between the heat transfer plates 
was fixed at 10°C and the plate separation distance at 
L = 40 mm (Ra = 5 x 10 5 ). The cavity aspect ratio was 5 
in both the horizontal and vertical planes (H/L and B/L). 
Morrison and Tran measured the velocity components by 
using Laser-Doppler anemometry. The nondimensional 
z-component velocity is compared between the present 
results and Morrison and Tran’s experimental data in 
figures 6 and 7. Good agreement is observed, indicating 
that the Boussinesq approximation is adequate for this 
type of flow. 

Concluding Remarks 

In the present study, it is shown that the Boussinesq 
approximation is valid for the analysis of heat transfer 
in an incompressible medium with mild temperature 
variations where radiation or boiling heat transfer can be 
neglected. The resulting formulation is validated using a 
version of the INS3D code. Computed results show that 
forced, mixed, and natural convection problems can be 
accurately predicted using the Boussinesq approximation, 
even at Ra = 5 x 10 5 in the case of Morrison’s experiment 
(ref. 15). This indicates that the interaction between 
pseudocompressibility and buoyancy force is properly 
accounted for by the present method. Overall, the cost of 
computing attributed to the temperature equation has been 
increased less than 5 percent. 
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temperature stream function vorticity 


(a) Forced convection (Re = 20, Gr = 0) 



temperature stream function vorticity 


(b) Mixed convection (Re = 20, Gr = 1600) 



temperature stream function vorticity 


(c) Natural convection (Ra = 10 3 ) 

Figure 2. Flow around a heated circular cylinder: Comparison of incompressible Navier-Stokes (INS) and stream function- 
vorticity (SV) computations. 
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(a) Forced convection 



(b) Mixed convection 



(c) Natural convection 


Figure 3. Heat transfer coefficient for flow around a heated circular cylinder: Comparison of computed results and 
experimental data. 
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(b) Aspect ratio of 2.0 (Ra = 10 5 : type A of bifurcation flow) 



INS SV 


INS SV 




(c) Aspect ratio of 2.0 (Ra = 10 5 : type B of bifurcation flow) 


Figure 4. Thermally driven bifurcation in two-dimensional rectangular cavity: Contour plot of the temperature, the stream 
function, and the vorticity. 
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Figure 7. Vertical velocity component, w, at y/H = 0.5 for natural convection problem in a 3-D cavity: Tp~T c = 10°C, 
Ra = 10 5 , z* = z/H. 
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